Binary (Logistic) Regression with Two Quantitative Predictors

To illustrate a regression model with two quantitative predictors, we will develop a logistic regression model for obesity as a function of age (adults 18-79) and household income (up to $100,000) for the US population. We will treat household income, an ordinal variable, as a quantitative variable coded using the midpoint values for each category.

Age and household income are quantitative variables. As we’ve discussed, we typically model the effects of quantitative variables using natural cubic splines. In formulating our initial model, we need to specify (1) the number of parameters (df) (complexity or degree of nonlinearity) for the age effect, (2) the number of parameters (df) (complexity or degree of nonlinearity) for the income effect, and (3) whether to assume additivity or to allow for possible interaction between age and income.

As we’ve discussed, a natural cubic spline with 4 df typically represents an appropriate balance between flexibility and loss of precision caused by overfitting. Small samples (\(N < 30\)) may require the use of 2 df, and moderate samples (\(N < 100\)) may require the use of 3 df. Empirical evidence suggests that more than 4 df are seldom required in a natural cubic spline regression model. Cubic spline models with 5 or 6 df would typically only be considered if the sample size was extremely large.

Based on these considerations, we will represent age and income using a natural cubic spline with 4 df.

As discussed previously for two categorical predictors, the choice of an additive model or an interaction model depends primarily on the scientific context and, to a lesser extent, on the available sample size. In smaller samples, the number of parameters \(p\) (regression coefficients) that can be reliably estimated is \(m/15\), where \(m\) is the effective sample size (for logistic regression, \(m\) is the minimum of the number of events and the number of non-events in our sample).

In this specific case, in my opinion, an additive model for the effects of age and income is potentially reasonable.

So, our initial model is an additive logistic regression model for prevalence of obesity as a function of age (represented as a natural cubic spline with 4 df) and income (represented as a natural cubic spline with 4 df).

Additive Model

We denote the obesity outcome for subject \(i\) by \(y_i\), the age for subject \(i\) by \(x_{1i}\), and the income for subject \(i\) by \(x_{2i}\).

An additive logistic regression model for (prevalence of) obesity as a function of age and income is \[ y_i \sim \mbox{Bernoulli}\left\{p_i\right\} \\ \mbox{logit}\left\{p_i\right\} = \log{\left\{\frac{p_i}{1-p_i}\right\}} = \alpha_0 + f_{a}(x_{1i},{\bf \beta}) + f_{i}(x_{2i},{\bf \gamma})\]

For modeling, age will be represented by a basis for the natural cubic spline with 4 df (5 knots), e.g. \(a_1(x_{2}), a_2(x_{2}), a_3(x_{2}), a_4(x_{2})\).

In our regression model, \[f_{a}(x_{2},{\bf \gamma}) = \gamma_1 a_1(x_{2}) + \gamma_2 a_2(x_{2}) + \gamma_3 a_3(x_{2}) + \gamma_4 a_4(x_{2})\]

For modeling, income will be represented by a basis for the natural cubic spline with 4 df (5 knots), e.g. \(i_1(x_{2}), i_2(x_{2}), i_3(x_{2}), i_4(x_{2})\).

In our regression model, \[f_{i}(x_{2},{\bf \gamma}) = \gamma_1 i_1(x_{2}) + \gamma_2 i_2(x_{2}) + \gamma_3 i_3(x_{2}) + \gamma_4 i_4(x_{2})\]

LIBNAME NHANES "../../../BMI552_Datasets/NHANES_May2023/";
*LIBNAME NHANES "~/my_shared_file_links/u48718377/NHANES";
OPTIONS FMTSEARCH=(NHANES.formats_NHANES work) NODATE NONUMBER;
TITLE1; TITLE2;
ODS NOPROCTITLE;
ODS GRAPHICS / LOESSMAXOBS=20000;

DATA Ron;
  SET NHANES.NHANES;
RUN;

ODS SELECT NONE;
PROC LOGISTIC DATA=Ron;
  ODS SELECT ResponseProfile;
*  ODS SELECT SplineKnots;
  EFFECT AgeCS=SPLINE(Age / NATURALCUBIC BASIS=TPF(NOINT) DEGREE=3 
                            KNOTMETHOD=PERCENTILELIST(5 27.5 50 72.5 95) DETAILS);
  EFFECT IncomeCS=SPLINE(HHIncomeMid / NATURALCUBIC BASIS=TPF(NOINT) DEGREE=3 
                            KNOTMETHOD=PERCENTILELIST(5 27.5 50 72.5 95) DETAILS);
    
  MODEL Obese(EVENT="Yes") = AgeCS IncomeCS / LINK=LOGIT;

  OUTPUT OUT=Tmp RESCHI=Residual PREDICTED=Fitted; 

  STORE AgeIncome;
RUN;
ODS SELECT ALL;

PROC FORMAT;
  VALUE IncomeGrp
    LOW-19999="0-19999"
    20000-39999="20000-39999"
    40000-59999="40000-59999"
    60000-79999="60000-79999"
    80000-99999="80000-99999"
    ;

  VALUE AgeGrp
    LOW-29="18-29"
    30-39="30-39"
    40-49="40-49"
    50-59="50-59"
    60-69="60-69"
    70-79="70-79"
    ;
RUN;
        
PROC SGPANEL DATA=Tmp;
  FORMAT HHIncomeMid IncomeGrp.;
  LABEL HHIncomeMid="Household Income";
  PANELBY HHIncomeMid / LAYOUT=ROWLATTICE;
    REFLINE 0 / AXIS=Y;
  LOESS X=Age Y=Residual / JITTER; 
RUN;

PROC SGPANEL DATA=Tmp;
  FORMAT Age AgeGrp.;
  PANELBY Age / LAYOUT=ROWLATTICE;
    REFLINE 0 / AXIS=Y;
  LOESS X=HHIncomeMid Y=Residual / JITTER; 
RUN;